Advanced visualization DZ

Tasks 1 and 2

1. Загрузите датасет very_low_birthweight.RDS (лежит в папке домашнего задания). 

Сделайте копию датасета, в которой удалите колонки с количеством пропусков больше 100, а затем удалите все строки с пропусками. 

2. Постройте графики плотности распределения для числовых переменных. Удалите выбросы, если таковые имеются. Преобразуйте категориальные переменные в факторы. Для любых двух числовых переменных раскрасьте график по переменной ‘inout’.

ds <- readRDS('very_low_birthweight.RDS')
glimpse(ds)

## Rows: 671
## Columns: 26
## $ birth    <dbl> 81.511, 81.514, 81.552, 81.558, 81.593, 81.602, 81.610, NA, N…
## $ exit     <dbl> 81.604, 81.539, 81.552, 81.667, 81.599, 81.771, 81.697, NA, N…
## $ hospstay <int> 34, 9, -2, 40, 2, 62, 32, NA, NA, 28, 38, NA, 62, 69, 1, 93, …
## $ lowph    <dbl> NA, 7.250000, 7.059998, 7.250000, 6.969997, 7.189999, 7.32000…
## $ pltct    <int> 100, 244, 114, 182, 54, NA, 282, NA, NA, 153, 229, NA, 182, 3…
## $ race     <fct> white, white, black, black, black, white, black, NA, NA, blac…
## $ bwt      <int> 1250, 1370, 620, 1480, 925, 940, 1255, 600, 700, 1350, 1310, …
## $ gest     <dbl> 35.0, 32.0, 23.0, 32.0, 28.0, 28.0, 29.5, 26.0, 24.0, 34.0, 3…
## $ inout    <fct> born at Duke, born at Duke, born at Duke, born at Duke, born …
## $ twn      <int> 0, 0, 0, 0, 0, 0, 0, NA, NA, 0, 0, NA, 0, 0, 0, 0, 0, 1, 1, 0…
## $ lol      <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ magsulf  <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ meth     <int> 0, 1, 0, 1, 0, 1, 1, NA, NA, 1, 0, NA, 0, 0, 1, 1, 1, 1, 1, 1…
## $ toc      <int> 0, 0, 1, 0, 0, 0, 0, NA, NA, 0, 0, NA, 0, 0, 1, 1, 0, 0, 0, 1…
## $ delivery <fct> abdominal, abdominal, vaginal, vaginal, abdominal, abdominal,…
## $ apg1     <int> 8, 7, 1, 8, 5, 8, 9, NA, NA, 4, 6, NA, 6, 6, 2, 4, 8, 7, 1, 8…
## $ vent     <int> 0, 0, 1, 0, 1, 1, 0, NA, NA, 0, 1, NA, 0, 0, 1, 1, 0, 0, 1, 1…
## $ pneumo   <int> 0, 0, 0, 0, 1, 0, 0, NA, NA, 0, 0, NA, 1, 0, 1, 0, 0, 0, 1, 1…
## $ pda      <int> 0, 0, 0, 0, 0, 0, 0, NA, NA, 0, 0, NA, 0, 0, 0, 0, 0, 0, 0, 0…
## $ cld      <int> 0, 0, NA, 0, 0, 0, 0, NA, NA, 0, 0, NA, 1, 0, 0, 1, 0, 0, 0, …
## $ pvh      <fct> NA, NA, NA, NA, definite, absent, NA, NA, absent, NA, NA, NA,…
## $ ivh      <fct> NA, NA, NA, NA, definite, absent, NA, NA, absent, NA, NA, NA,…
## $ ipe      <fct> NA, NA, NA, NA, NA, absent, NA, NA, absent, NA, NA, NA, absen…
## $ year     <dbl> 81.51196, 81.51471, 81.55304, 81.55847, 81.59406, 81.60229, 8…
## $ sex      <fct> female, female, female, male, female, female, female, NA, NA,…
## $ dead     <int> 0, 0, 1, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0…

summary(ds)

##      birth            exit          hospstay            lowph      
##  Min.   :81.51   Min.   :68.53   Min.   :-6574.00   Min.   :6.530  
##  1st Qu.:83.52   1st Qu.:83.58   1st Qu.:   16.00   1st Qu.:7.130  
##  Median :84.90   Median :84.96   Median :   37.00   Median :7.210  
##  Mean   :84.75   Mean   :84.84   Mean   :   40.36   Mean   :7.202  
##  3rd Qu.:86.07   3rd Qu.:86.17   3rd Qu.:   62.00   3rd Qu.:7.310  
##  Max.   :87.48   Max.   :96.87   Max.   : 3668.00   Max.   :7.550  
##  NA's   :21      NA's   :31      NA's   :31         NA's   :62     
##      pltct                    race          bwt            gest      
##  Min.   : 16.0   black          :369   Min.   : 400   Min.   :22.00  
##  1st Qu.:143.0   native American: 16   1st Qu.: 900   1st Qu.:27.00  
##  Median :202.0   oriental       :  4   Median :1120   Median :29.00  
##  Mean   :201.6   white          :257   Mean   :1094   Mean   :28.87  
##  3rd Qu.:252.0   NA's           : 25   3rd Qu.:1310   3rd Qu.:31.00  
##  Max.   :571.0                         Max.   :1580   Max.   :40.00  
##  NA's   :70                            NA's   :2      NA's   :4      
##           inout          twn              lol             magsulf      
##  born at Duke:547   Min.   :0.0000   Min.   :  0.000   Min.   :0.0000  
##  transported :121   1st Qu.:0.0000   1st Qu.:  0.000   1st Qu.:0.0000  
##  NA's        :  3   Median :0.0000   Median :  3.500   Median :0.0000  
##                     Mean   :0.2074   Mean   :  8.438   Mean   :0.1344  
##                     3rd Qu.:0.0000   3rd Qu.:  9.000   3rd Qu.:0.0000  
##                     Max.   :1.0000   Max.   :192.000   Max.   :1.0000  
##                     NA's   :20       NA's   :381       NA's   :247     
##       meth             toc              delivery        apg1      
##  Min.   :0.0000   Min.   :0.0000   abdominal:314   Min.   :0.000  
##  1st Qu.:0.0000   1st Qu.:0.0000   vaginal  :335   1st Qu.:2.000  
##  Median :0.0000   Median :0.0000   NA's     : 22   Median :5.000  
##  Mean   :0.4372   Mean   :0.2248                   Mean   :4.903  
##  3rd Qu.:1.0000   3rd Qu.:0.0000                   3rd Qu.:7.000  
##  Max.   :1.0000   Max.   :1.0000                   Max.   :9.000  
##  NA's   :106      NA's   :106                      NA's   :34     
##       vent            pneumo            pda              cld        
##  Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   Min.   :0.0000  
##  1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.0000  
##  Median :1.0000   Median :0.0000   Median :0.0000   Median :0.0000  
##  Mean   :0.5803   Mean   :0.1969   Mean   :0.2087   Mean   :0.2694  
##  3rd Qu.:1.0000   3rd Qu.:0.0000   3rd Qu.:0.0000   3rd Qu.:1.0000  
##  Max.   :1.0000   Max.   :1.0000   Max.   :1.0000   Max.   :1.0000  
##  NA's   :30       NA's   :26       NA's   :29       NA's   :66      
##        pvh            ivh            ipe           year           sex     
##  absent  :360   absent  :442   absent  :472   Min.   :81.51   female:320  
##  definite:125   definite: 75   definite: 38   1st Qu.:83.52   male  :330  
##  possible: 41   possible: 10   possible: 17   Median :84.91   NA's  : 21  
##  NA's    :145   NA's    :144   NA's    :144   Mean   :84.76               
##                                               3rd Qu.:86.07               
##                                               Max.   :87.48               
##                                               NA's   :21                  
##       dead       
##  Min.   :0.0000  
##  1st Qu.:0.0000  
##  Median :0.0000  
##  Mean   :0.2146  
##  3rd Qu.:0.0000  
##  Max.   :1.0000  
## 

factors <- c("race","inout","delivery", "sex", "twn", "vent", "pneumo", "pda", "cld", "dead")

ds_filtered0 <-  ds %>%
  #удаление колонок где больше 100 пропусков
  dplyr::select(where(~sum(is.na(.x)) <= 100)) %>%
  
  #перевод категориальных переменных в факторы
  mutate(across(all_of(factors), as.factor)) %>%
  
  #перекодировка факторов
  mutate(vent = recode(vent, '0' = "No ventilation", '1' = "Ventilation"),
         pneumo = recode(pneumo, '0' = 'No pneumothorax', '1' = 'Pneumothorax'),
         pda = recode(pda, '0' = 'No patent ductus arteriosus', '1' = 'Patent ductus arteriosus'),
         cld = recode(cld, '0' = 'No oxygen', '1' = 'On oxygen'),
         dead = recode(dead, '0' = 'Alive', '1' = 'Dead'),
         ID =as.factor(row_number())) %>%
  
  #удаление выбросов
  filter(if_any(where(is.numeric), 
                ~ .x >= quantile(.x, 0.25, na.rm = TRUE) - 1.5 * IQR(.x, na.rm = TRUE) & 
                  .x <= quantile(.x, 0.75, na.rm = TRUE) + 1.5 * IQR(.x, na.rm = TRUE))) 

#удаление стоблцов с NA
ds_filtered <- ds_filtered0 %>% drop_na
  

summary(ds_filtered)

##      birth            exit          hospstay           lowph      
##  Min.   :81.51   Min.   :81.05   Min.   :-295.00   Min.   :6.530  
##  1st Qu.:83.43   1st Qu.:83.56   1st Qu.:  21.00   1st Qu.:7.135  
##  Median :84.77   Median :84.87   Median :  40.00   Median :7.220  
##  Mean   :84.63   Mean   :84.76   Mean   :  47.04   Mean   :7.215  
##  3rd Qu.:85.83   3rd Qu.:85.99   3rd Qu.:  64.00   3rd Qu.:7.320  
##  Max.   :87.48   Max.   :87.72   Max.   : 797.00   Max.   :7.550  
##                                                                   
##      pltct                    race          bwt            gest      
##  Min.   : 16.0   black          :303   Min.   : 400   Min.   :23.00  
##  1st Qu.:148.0   native American: 13   1st Qu.: 960   1st Qu.:28.00  
##  Median :204.0   oriental       :  4   Median :1160   Median :29.00  
##  Mean   :204.5   white          :211   Mean   :1136   Mean   :29.25  
##  3rd Qu.:256.0                         3rd Qu.:1330   3rd Qu.:31.00  
##  Max.   :571.0                         Max.   :1500   Max.   :36.00  
##                                                                      
##           inout     twn          delivery        apg1      
##  born at Duke:448   0:422   abdominal:262   Min.   :0.000  
##  transported : 83   1:109   vaginal  :269   1st Qu.:2.000  
##                                             Median :5.000  
##                                             Mean   :5.024  
##                                             3rd Qu.:7.000  
##                                             Max.   :9.000  
##                                                            
##              vent                 pneumo                             pda     
##  No ventilation:243   No pneumothorax:438   No patent ductus arteriosus:425  
##  Ventilation   :288   Pneumothorax   : 93   Patent ductus arteriosus   :106  
##                                                                              
##                                                                              
##                                                                              
##                                                                              
##                                                                              
##         cld           year           sex         dead           ID     
##  No oxygen:393   Min.   :81.51   female:264   Alive:467   2      :  1  
##  On oxygen:138   1st Qu.:83.43   male  :267   Dead : 64   4      :  1  
##                  Median :84.77                            5      :  1  
##                  Mean   :84.63                            7      :  1  
##                  3rd Qu.:85.83                            10     :  1  
##                  Max.   :87.48                            11     :  1  
##                                                           (Other):525

ds_filtered %>%
  pivot_longer(cols = where(is.numeric), names_to = "variable", values_to = "value") %>%
  ggplot(aes(x = value, fill = inout)) +
  geom_density(alpha = 0.5) +
  facet_wrap(~ variable, scales = "free") +
  labs(title = "Density Plots of Numeric Variables", x = "Value", y = "Density") +
  theme_minimal()

Task 3

Проведите тест на сравнение значений колонки ‘lowph’ между группами в переменной inout. Вид статистического теста определите самостоятельно. Визуализируйте результат через библиотеку ‘rstatix’. Как бы вы интерпретировали результат, если бы знали, что более низкое значение lowph ассоциировано с более низкой выживаемостью?

t_test_result <- ds_filtered %>%
  t_test(lowph ~ inout) %>%
  add_significance() %>%
  mutate(y.position = max(ds_filtered$lowph, na.rm = TRUE) * 1.05) 

t_test_result

## # A tibble: 1 × 10
##   .y.   group1    group2    n1    n2 statistic    df       p p.signif y.position
##   <chr> <chr>     <chr>  <int> <int>     <dbl> <dbl>   <dbl> <chr>         <dbl>
## 1 lowph born at … trans…   448    83      5.57  111. 1.77e-7 ****           7.93

ds_filtered %>%
  ggplot(aes(x = inout, y = lowph)) +
  geom_boxplot() +
  stat_pvalue_manual(t_test_result, label = "p") +
  labs(x = "Inout Group",
       y = "Low pH") +
  theme_minimal()

Я выбрала т-тест поскольку довольно много наблюдений, его результаты показывают, что уровень pH значимо отличается у детей родившихся в Дьюке и привезенных из других госпиталей. Это может говорить о том, что транспортировка детей ухудшает их прогноз, возможно им не успевают оказать необходимую помощь.

Task 4

Сделайте новый датафрейм, в котором оставьте только континуальные или ранговые данные, кроме ‘birth’, ‘year’ и ‘exit’. Сделайте корреляционный анализ этих данных. Постройте два любых типа графиков для визуализации корреляций.

continuous_ranked_vars <- ds_filtered %>%
  dplyr::select(where(~ is.numeric(.x) || is.ordered(.x))) %>%
  dplyr::select(-birth, -year, -exit)

glimpse(continuous_ranked_vars)

## Rows: 531
## Columns: 6
## $ hospstay <int> 9, 40, 2, 32, 28, 38, 62, 69, 1, 93, 44, 66, 65, 44, 70, 85, …
## $ lowph    <dbl> 7.250000, 7.250000, 6.969997, 7.320000, 7.160000, 7.039997, 7…
## $ pltct    <int> 244, 182, 54, 282, 153, 229, 182, 361, 378, 255, 186, 260, 18…
## $ bwt      <int> 1370, 1480, 925, 1255, 1350, 1310, 1110, 1180, 970, 770, 1490…
## $ gest     <dbl> 32.0, 32.0, 28.0, 29.5, 34.0, 32.0, 28.0, 28.0, 28.0, 26.0, 3…
## $ apg1     <int> 7, 8, 5, 9, 4, 6, 6, 6, 2, 4, 8, 1, 8, 5, 9, 9, 9, 6, 2, 1, 6…

cor_matrix <- cor(continuous_ranked_vars, use = "pairwise.complete.obs", method = "spearman")
cor_matrix

##            hospstay      lowph       pltct        bwt        gest       apg1
## hospstay  1.0000000 -0.2337033 -0.15825443 -0.4447105 -0.38014775 -0.1255054
## lowph    -0.2337033  1.0000000  0.29144157  0.3234317  0.36957021  0.2730238
## pltct    -0.1582544  0.2914416  1.00000000  0.2465972  0.08007581  0.2982137
## bwt      -0.4447105  0.3234317  0.24659722  1.0000000  0.71038991  0.3182650
## gest     -0.3801477  0.3695702  0.08007581  0.7103899  1.00000000  0.2674875
## apg1     -0.1255054  0.2730238  0.29821365  0.3182650  0.26748752  1.0000000

network_plot(cor_matrix)

chart.Correlation(continuous_ranked_vars, histogram = TRUE, pch = 19)

## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter

Task 5

Постройте иерархическую кластеризацию на этом датафрейме.

distance_matrix <- dist(continuous_ranked_vars, method = "euclidean")
hclust_result <- hclust(distance_matrix, method = "ward.D2")

label_colors <- ifelse(ds_filtered$dead == "Alive", "blue", "red")
fviz_dend(hclust_result, 
          show_labels = TRUE,
          label_cols = label_colors,
          main = "Hierarchical Clustering Dendrogram",
          sub = "", xlab = "", ylab = '')

## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## ℹ The deprecated feature was likely used in the factoextra package.
##   Please report the issue at <https://github.com/kassambara/factoextra/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Task 6

Сделайте одновременный график heatmap и иерархической кластеризации. Интерпретируйте результат.

pheatmap(cor_matrix, 
         clustering_distance_rows = "euclidean", 
         clustering_distance_cols = "euclidean", 
         clustering_method = "ward.D2", 
         display_numbers = TRUE, 
         number_format = "%.2f",
         main = "Heatmap with Hierarchical Clustering")

По данному графику видно, что числовые переменные можно разделить на 2 группы, которые между собой связаны внутри групп, не считая времени проведенного в госпитале, которая обратно коррелирует со всеми остальными переменными. В первой группе вес при рождении и продолжительность беременности, во второй - pH, значение по шкале Апгар и количество тромбоцитов.

Tasks 7-10

Проведите PCA анализ на этих данных. Проинтерпретируйте результат. Нужно ли применять шкалирование для этих данных перед проведением PCA?

Постройте biplot график для PCA. Раскрасьте его по значению колонки ‘dead’.

Переведите последний график в ‘plotly’. При наведении на точку нужно, чтобы отображалось id пациента.

data_scaled <- scale(continuous_ranked_vars)
pca_result <- PCA(data_scaled, graph = FALSE)

В данном случае данные нужно обязательно шкалировать, так как переменные в разных единицах.

Biplot

pca_ind <- as.data.frame(pca_result$ind$coord) %>% 
  mutate(dead = ds_filtered$dead, id = as.character(ds_filtered$ID)) 

pca_var <- as.data.frame(pca_result$var$coord)
pca_var$Variable <- rownames(pca_var)

p <- ggplot() +
  geom_point(data = pca_ind, aes(x = Dim.1, y = Dim.2, color = dead, text = id)) +
  geom_segment(data = pca_var, aes(x = 0, y = 0, xend = Dim.1, yend = Dim.2)) +
  geom_text(data = pca_var, aes(x = Dim.1, y = Dim.2, label = Variable), hjust = 1.2, vjust = 1.2) +
  scale_color_manual(values = c("blue", "red")) +
  labs(title = "Interactive PCA Biplot", x = "PC1", y = "PC2") +
  theme_minimal()

## Warning in geom_point(data = pca_ind, aes(x = Dim.1, y = Dim.2, color = dead, :
## Ignoring unknown aesthetics: text

ggplotly(p, tooltip = "text")

сложно интерпретировать этот график, я пыталась по разному менять размеры стрелок и точек, но для этих данных проще построить отдельный график со стрелками чтобы посмотреть на вклад переменных в главные компоненты

fviz_pca_var(pca_result, col.var = "contrib", gradient.cols = c("blue", "yellow", "red"),
             repel = TRUE, title = "Contribution of Variables to PCA")

На этом графике еще более четко видно разделение переменных по 3 группам описанным выше. Видно, что есть 2 переменные - вес и срок беременности, которые дают самый весомый вклад в первую компоненту и она объясняет довольно много дисперсии (почти 40%). Вероятно, это 2 главных фактора, которые влияют на состояние недоношенных младенцев. Другие 3 переменные, которые харакетризиют жизненные показатели младенцев тоже вносят большой вклад в первую компоненту, возможно они свзяны с осложнениями при беременности и родах, таких как преэклапсия или инфекции. Вторая компонента основана главным образом на переменной время проведенное в госпитале, вторая компонента объясняет мало дисперсии данных по сравнению с первой. Возможно потому что время пребывания не всегда связано с состоянием пациента, или эта связь сложная, к примеру мало прибывают те, кто умирает или наоборот в хорошем состоянии, а те, кому нужна помощь будут находиться дольше. Похоже на то, что первая компонента разделяет младенцев на умерших и выживших.

Некорректно использовать колонку dead, так как она не дает информации о том когда именно была смерть.

explained_variance <- pca_result$eig[, 2]
explained_variance_df <- data.frame(
  Component = paste0("PC", 1:length(explained_variance)),
  Variance = explained_variance
)

ggplot(explained_variance_df, aes(x = Component, y = Variance)) +
  geom_bar(stat = "identity", fill = "steelblue") +
  geom_text(aes(label = paste0(round(Variance, 1), "%")), vjust = -0.5, size = 4) +
  labs(
    title = "Explained Variance by Principal Components",
    x = "Principal Component",
    y = "Explained Variance (%)"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Еще хороший график чтобы понять какие кроме первых двух компонент могут еще вносить вклад. В нашем случае не супер информативно, тк мало переменных в исходном датасете.

Task 11

umap_result <- umap(continuous_ranked_vars, n_neighbors = 15, min_dist = 0.1, n_components = 2)
umap_data <- as.data.frame(umap_result$layout)
umap_data$dead <- ds_filtered$dead

ggplot(umap_data, aes(x = V1, y = V2, color = dead)) +
  geom_point(size = 3, alpha = 0.7) +
  labs(title = "UMAP Projection", x = "UMAP 1", y = "UMAP 2") +
  theme_minimal()

С помощью UMAP данные лучше разделились по 2 компоненте, чем по первой, умершие попали так же преимущественно в один кластер как и на PCA.

Task 12

Чтобы проверить как влияют параметры я сделала 4 варианта параметров и запустила для них UMAP

umap_configs <- list(
  list(n_neighbors = 5, min_dist = 0.1),
  list(n_neighbors = 15, min_dist = 0.1),
  list(n_neighbors = 15, min_dist = 0.5),
  list(n_neighbors = 30, min_dist = 0.5)
)

apply_umap <- function(config, data) {
  result <- umap(data, n_neighbors = config$n_neighbors, min_dist = config$min_dist, n_components = 2)
  as.data.frame(result$layout) %>%
    mutate(n_neighbors = config$n_neighbors, min_dist = config$min_dist)
}

umap_results <- bind_rows(lapply(umap_configs, apply_umap, data = continuous_ranked_vars))

umap_results <- umap_results %>%
  mutate(dead = rep(ds_filtered$dead, times = length(umap_configs)))


ggplot(umap_results, aes(x = V1, y = V2, color = dead)) +
  geom_point(size = 2, alpha = 0.7) +
  facet_wrap(~ n_neighbors + min_dist, nrow = 2, ncol = 2, labeller = label_both) +
  labs(title = "UMAP Projections with Different Parameters",
       x = "UMAP 1", y = "UMAP 2") +
  scale_color_manual(values = c("blue", "red")) +
  theme_minimal()

Проекции с разными параметрами имеют разное количество кластеров: чем больше оба параметра - тем меньше кластеров и тем лучше видно как точки внутри кластеров расположены относительно друг друга. На самом первом графике максимальное количество групп точек, так как алгоритм учитывает только 5 соседних наблюдений, на самом же последнем графике точки все точки сгруппированы вместе.

Task 13

Пермутация

set.seed(123)
data_perm_50 <- continuous_ranked_vars %>%
  mutate(bwt = ifelse(runif(n()) <= 0.5, sample(bwt), bwt))

data_perm_100 <- continuous_ranked_vars %>%
  mutate(bwt = sample(bwt))

Функции для PCA и UMAP

run_pca <- function(data) {
  pca_result <- PCA(data, graph = FALSE)
  return(pca_result)
}

run_umap <- function(data) {
  umap_result <- umap(data, n_neighbors = 15, min_dist = 0.1, n_components = 2)
  return(as.data.frame(umap_result$layout))
}

#PCA и UMAP для оригинальных данных
pca_original <- run_pca(continuous_ranked_vars)
umap_original <- run_umap(continuous_ranked_vars)

#PCA и UMAP для данных с пермутацией 50%
pca_perm_50 <- run_pca(data_perm_50)
umap_perm_50 <- run_umap(data_perm_50)

#PCA и UMAP для данных с пермутацией 100%
pca_perm_100 <- run_pca(data_perm_100)
umap_perm_100 <- run_umap(data_perm_100)

explained_variance <- tibble(
  Component = paste0("PC", 1:length(pca_original$eig[, 2])),
  Original = pca_original$eig[, 2],
  Permuted_50 = pca_perm_50$eig[, 2] ,
  Permuted_100 = pca_perm_100$eig[, 2]
)

print(explained_variance)

## # A tibble: 6 × 4
##   Component Original Permuted_50 Permuted_100
##   <chr>        <dbl>       <dbl>        <dbl>
## 1 PC1          39.4        32.7         30.3 
## 2 PC2          17.9        17.4         18.0 
## 3 PC3          15.0        15.3         16.3 
## 4 PC4          12.2        14.1         14.3 
## 5 PC5          11.0        12.0         12.0 
## 6 PC6           4.47        8.45         9.03

Из таблицы выше видно, что пермутация данных ухудшила PCA - первая компонента стала объяснять меньше вариабельности почти на 25% для 100% пермутированных данных по сравнению с оригинальными.

pca_data <- list(
  Original = as.data.frame(pca_original$ind$coord) %>% mutate(Method = "PCA", Data = "Original", dead = ds_filtered$dead),
  Permuted_50 = as.data.frame(pca_perm_50$ind$coord) %>% mutate(Method = "PCA", Data = "Permuted 50%", dead = ds_filtered$dead),
  Permuted_100 = as.data.frame(pca_perm_100$ind$coord) %>% mutate(Method = "PCA", Data = "Permuted 100%", dead = ds_filtered$dead)
) %>% bind_rows()

umap_data <- list(
  Original = umap_original %>% 
    mutate(Method = "UMAP", Data = "Original", dead = ds_filtered$dead),
  
  Permuted_50 = umap_perm_50 %>% 
    mutate(Method = "UMAP", Data = "Permuted 50%", dead = ds_filtered$dead),
  
  Permuted_100 = umap_perm_100 %>% 
    mutate(Method = "UMAP", Data = "Permuted 100%", dead = ds_filtered$dead)
) %>% bind_rows()


combined_data <- bind_rows(
  pca_data %>% rename(Dim1 = Dim.1, Dim2 = Dim.2),
  umap_data %>% rename(Dim1 = V1, Dim2 = V2)
)

ggplot(combined_data, aes(x = Dim1, y = Dim2, color = dead)) +
  geom_point(size = 2, alpha = 0.5) +
  facet_grid(Method ~ Data) + 
  labs(title = "PCA and UMAP Projections with Permuted Data",
       x = "Dimension 1", y = "Dimension 2") +
  theme_minimal() +
  theme(legend.position = "bottom")

Честно говоря на этих данных сложно увидеть эффект пермутации, в теории кластеры должны были значительно поменяться, но на PCA эффект не очень заметен, на UMAP проекции с пермутациями отличаются от изначальной тем, что точки для мертвых не группируются вместе. Чтобы лучше увидеть эффект для PCA нарисовала еще графики со стрелками, на них видно, что bwt вносит меньше вклада чем больше пермутированных значений.

pca_var_original <- fviz_pca_var(pca_original, 
                                 col.var = "contrib", 
                                 gradient.cols = c("blue", "green", "red"),
                                 repel = TRUE, 
                                 title = "PCA Variables: Original Data")

# Permuted 50% Data
pca_var_perm_50 <- fviz_pca_var(pca_perm_50, 
                                col.var = "contrib", 
                                gradient.cols = c("blue", "green", "red"),
                                repel = TRUE, 
                                title = "PCA Variables: Permuted 50%")

# Permuted 100% Data
pca_var_perm_100 <- fviz_pca_var(pca_perm_100, 
                                 col.var = "contrib", 
                                 gradient.cols = c("blue", "green", "red"),
                                 repel = TRUE, 
                                 title = "PCA Variables: Permuted 100%")


grid.arrange(pca_var_original, pca_var_perm_50, pca_var_perm_100, nrow = 1)

Task 14

Давайте проведем анализ чувствительности. Проведите анализ, как в шагах 4-6 для оригинального с удалением всех строк с пустыми значениями (т.е. включая колонки с количеством пропущенных значений больше 100), а затем для оригинального датафрейма с импутированием пустых значений средним или медианой. Как отличаются получившиеся результаты? В чем преимущества и недостатки каждого подхода?

#импутировала NA регрессией с учетом других значений
vars_to_exclude <- c("ID", "year", "birth", "exit")

ds_filtered_imp0 <- ds_filtered0 %>% 
  dplyr::select(-all_of(vars_to_exclude)) %>%
  mice() %>%
  complete()

## 
##  iter imp variable
##   1   1  hospstay  lowph  pltct  race  bwt  gest  inout  twn  delivery  apg1  vent  pneumo  pda  cld  sex
##   1   2  hospstay  lowph  pltct  race  bwt  gest  inout  twn  delivery  apg1  vent  pneumo  pda  cld  sex
##   1   3  hospstay  lowph  pltct  race  bwt  gest  inout  twn  delivery  apg1  vent  pneumo  pda  cld  sex
##   1   4  hospstay  lowph  pltct  race  bwt  gest  inout  twn  delivery  apg1  vent  pneumo  pda  cld  sex
##   1   5  hospstay  lowph  pltct  race  bwt  gest  inout  twn  delivery  apg1  vent  pneumo  pda  cld  sex
##   2   1  hospstay  lowph  pltct  race  bwt  gest  inout  twn  delivery  apg1  vent  pneumo  pda  cld  sex
##   2   2  hospstay  lowph  pltct  race  bwt  gest  inout  twn  delivery  apg1  vent  pneumo  pda  cld  sex
##   2   3  hospstay  lowph  pltct  race  bwt  gest  inout  twn  delivery  apg1  vent  pneumo  pda  cld  sex
##   2   4  hospstay  lowph  pltct  race  bwt  gest  inout  twn  delivery  apg1  vent  pneumo  pda  cld  sex
##   2   5  hospstay  lowph  pltct  race  bwt  gest  inout  twn  delivery  apg1  vent  pneumo  pda  cld  sex
##   3   1  hospstay  lowph  pltct  race  bwt  gest  inout  twn  delivery  apg1  vent  pneumo  pda  cld  sex
##   3   2  hospstay  lowph  pltct  race  bwt  gest  inout  twn  delivery  apg1  vent  pneumo  pda  cld  sex
##   3   3  hospstay  lowph  pltct  race  bwt  gest  inout  twn  delivery  apg1  vent  pneumo  pda  cld  sex
##   3   4  hospstay  lowph  pltct  race  bwt  gest  inout  twn  delivery  apg1  vent  pneumo  pda  cld  sex
##   3   5  hospstay  lowph  pltct  race  bwt  gest  inout  twn  delivery  apg1  vent  pneumo  pda  cld  sex
##   4   1  hospstay  lowph  pltct  race  bwt  gest  inout  twn  delivery  apg1  vent  pneumo  pda  cld  sex
##   4   2  hospstay  lowph  pltct  race  bwt  gest  inout  twn  delivery  apg1  vent  pneumo  pda  cld  sex
##   4   3  hospstay  lowph  pltct  race  bwt  gest  inout  twn  delivery  apg1  vent  pneumo  pda  cld  sex
##   4   4  hospstay  lowph  pltct  race  bwt  gest  inout  twn  delivery  apg1  vent  pneumo  pda  cld  sex
##   4   5  hospstay  lowph  pltct  race  bwt  gest  inout  twn  delivery  apg1  vent  pneumo  pda  cld  sex
##   5   1  hospstay  lowph  pltct  race  bwt  gest  inout  twn  delivery  apg1  vent  pneumo  pda  cld  sex
##   5   2  hospstay  lowph  pltct  race  bwt  gest  inout  twn  delivery  apg1  vent  pneumo  pda  cld  sex
##   5   3  hospstay  lowph  pltct  race  bwt  gest  inout  twn  delivery  apg1  vent  pneumo  pda  cld  sex
##   5   4  hospstay  lowph  pltct  race  bwt  gest  inout  twn  delivery  apg1  vent  pneumo  pda  cld  sex
##   5   5  hospstay  lowph  pltct  race  bwt  gest  inout  twn  delivery  apg1  vent  pneumo  pda  cld  sex

ds_filtered_imp <- bind_cols(ds_filtered_imp0, ds_filtered0 %>%  dplyr::select(all_of(vars_to_exclude)))

continuous_ranked_vars_imp <- ds_filtered_imp %>%
  dplyr::select(where(~ is.numeric(.x) || is.ordered(.x))) %>%
  dplyr::select(-birth, -year, -exit)

cor_matrix <- cor(continuous_ranked_vars_imp, use = "pairwise.complete.obs", method = "spearman")

network_plot(cor_matrix)

distance_matrix <- dist(continuous_ranked_vars_imp, method = "euclidean")
hclust_result <- hclust(distance_matrix, method = "ward.D2")
label_colors <- ifelse(ds_filtered0$dead == "Alive", "blue", "red")

fviz_dend(hclust_result, 
          show_labels = TRUE,
          label_cols = label_colors,
          main = "Hierarchical Clustering Dendrogram",
          sub = "", xlab = "", ylab = '')

pheatmap(cor_matrix, 
         clustering_distance_rows = "euclidean", 
         clustering_distance_cols = "euclidean", 
         clustering_method = "ward.D2", 
         display_numbers = TRUE, 
         number_format = "%.2f",
         main = "Heatmap with Hierarchical Clustering")

chart.Correlation(continuous_ranked_vars_imp, histogram = TRUE, pch = 19)

## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter

chart.Correlation(continuous_ranked_vars, histogram = TRUE, pch = 19)

## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter

Лучше всего видно как повлияла импутация на двух последних графиках: значения корреляции немного поменялись, для некоторых пар значимная корреляция есть только на одном из графиков. Наибольшие изменения для переменных gest и hospstay, у них даже форма гистограмм поменялась. Но в данном анализе сложно понять что больше повлияло - импутация или то что больше 100 строк данных добавилось.

Task 15

data_scaled <- scale(continuous_ranked_vars_imp)
pca_result <- PCA(data_scaled, graph = FALSE)

fviz_pca_var(pca_result, col.var = "contrib", gradient.cols = c("blue", "yellow", "red"),
             repel = TRUE, title = "Contribution of Variables to PCA")

pca_ind <- as.data.frame(pca_result$ind$coord) %>% 
  mutate(dead = ds_filtered0$dead, id = as.character(ds_filtered0$ID)) 

pca_var <- as.data.frame(pca_result$var$coord)
pca_var$Variable <- rownames(pca_var)

p <- ggplot() +
  geom_point(data = pca_ind, aes(x = Dim.1, y = Dim.2, color = dead, text = id)) +
  geom_segment(data = pca_var, aes(x = 0, y = 0, xend = Dim.1, yend = Dim.2)) +
  geom_text(data = pca_var, aes(x = Dim.1, y = Dim.2, label = Variable), hjust = 1.2, vjust = 1.2) +
  scale_color_manual(values = c("blue", "red")) +
  labs(title = "Interactive PCA Biplot", x = "PC1", y = "PC2") +
  theme_minimal()

## Warning in geom_point(data = pca_ind, aes(x = Dim.1, y = Dim.2, color = dead, :
## Ignoring unknown aesthetics: text

ggplotly(p, tooltip = "text")

umap_result <- umap(continuous_ranked_vars_imp, n_neighbors = 15, min_dist = 0.1, n_components = 2)
umap_data <- as.data.frame(umap_result$layout)
umap_data$dead <- ds_filtered0$dead

ggplot(umap_data, aes(x = V1, y = V2, color = dead)) +
  geom_point(size = 3, alpha = 0.7) +
  labs(title = "UMAP Projection", x = "UMAP 1", y = "UMAP 2") +
  theme_minimal()

Видно, что после импутации первая компонента PCA стала объяснять больше изменчивости и на графике проекции умершие и выжившие разделяются на группы гораздо лучше, что говорит о том, что по этим данным лучше получится определить связь смертности с указаными переменными. Еще интересно, что появились достаточно четкие выбросы на PCA, но их нет на UMAP. На UMAP по ощущению импутация вообще не повлияла.

Интересно, что в лекциях не было ничего сказано, что t-SNE, который похож на UMAP, но в биоинформатике гораздо более популярен. Попробовала им тоже сравнить данные до и после импутации и получилось довольно сильное отличие.

set.seed(42)
tsne_result <- Rtsne(as.matrix(continuous_ranked_vars), dims = 2, perplexity = 30, verbose = TRUE, max_iter = 500)

## Performing PCA
## Read the 531 x 6 data matrix successfully!
## OpenMP is working. 1 threads.
## Using no_dims = 2, perplexity = 30.000000, and theta = 0.500000
## Computing input similarities...
## Building tree...
## Done in 0.06 seconds (sparsity = 0.209277)!
## Learning embedding...
## Iteration 50: error is 56.340360 (50 iterations in 0.03 seconds)
## Iteration 100: error is 53.994127 (50 iterations in 0.03 seconds)
## Iteration 150: error is 53.890921 (50 iterations in 0.03 seconds)
## Iteration 200: error is 53.880156 (50 iterations in 0.03 seconds)
## Iteration 250: error is 53.876007 (50 iterations in 0.03 seconds)
## Iteration 300: error is 0.592579 (50 iterations in 0.03 seconds)
## Iteration 350: error is 0.507048 (50 iterations in 0.03 seconds)
## Iteration 400: error is 0.491431 (50 iterations in 0.03 seconds)
## Iteration 450: error is 0.482132 (50 iterations in 0.03 seconds)
## Iteration 500: error is 0.471446 (50 iterations in 0.03 seconds)
## Fitting performed in 0.29 seconds.

tsne_data <- data.frame(tsne_result$Y, dead = ds_filtered$dead)
colnames(tsne_data) <- c("Dim1", "Dim2", "dead")

# t-SNE визуализация
ggplot(tsne_data, aes(x = Dim1, y = Dim2, color = dead)) +
  geom_point(size = 2, alpha = 0.7) +
  labs(title = "t-SNE Visualization", x = "Dimension 1", y = "Dimension 2") +
  theme_minimal()

set.seed(42)
tsne_result <- Rtsne(as.matrix(continuous_ranked_vars_imp), dims = 2, perplexity = 30, verbose = TRUE, max_iter = 500)

## Performing PCA
## Read the 670 x 6 data matrix successfully!
## OpenMP is working. 1 threads.
## Using no_dims = 2, perplexity = 30.000000, and theta = 0.500000
## Computing input similarities...
## Building tree...
## Done in 0.09 seconds (sparsity = 0.166180)!
## Learning embedding...
## Iteration 50: error is 59.616393 (50 iterations in 0.06 seconds)
## Iteration 100: error is 55.926854 (50 iterations in 0.04 seconds)
## Iteration 150: error is 55.684807 (50 iterations in 0.04 seconds)
## Iteration 200: error is 55.658898 (50 iterations in 0.04 seconds)
## Iteration 250: error is 55.665700 (50 iterations in 0.04 seconds)
## Iteration 300: error is 0.714897 (50 iterations in 0.04 seconds)
## Iteration 350: error is 0.582656 (50 iterations in 0.04 seconds)
## Iteration 400: error is 0.546329 (50 iterations in 0.05 seconds)
## Iteration 450: error is 0.534019 (50 iterations in 0.05 seconds)
## Iteration 500: error is 0.523819 (50 iterations in 0.05 seconds)
## Fitting performed in 0.45 seconds.

tsne_data <- data.frame(tsne_result$Y, dead = ds_filtered0$dead)
colnames(tsne_data) <- c("Dim1", "Dim2", "dead")

# t-SNE визуализация
ggplot(tsne_data, aes(x = Dim1, y = Dim2, color = dead)) +
  geom_point(size = 2, alpha = 0.7) +
  labs(title = "t-SNE Visualization", x = "Dimension 1", y = "Dimension 2") +
  theme_minimal()